###
#
# Tyler Romualdi, Tyler Girard, Mathieu Turgeon, Yannick Dufresne, Takeshi Iida, and Tetsuya Matsubayashi. 
# "The Multidimensional Structure of Risk: How Dread and Control Shape Perceptions Toward Artificial Intelligence." 
# Journal of Risk Research. 
#
# This file imports the original survey data and prepares it for all subsequent analyses.
# The data for the predictive models is saved to the folder "2 Clean Data"
# The output figures are saved to the folder "4 Figures"
# The output tables are saved to the folder "5 Tables"
#
#
###

#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#
#### Maintain Folder Structure Using "Here" Package ####
#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#

library(here)
here()

#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#
#### Libraries ####
#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#

library(tidyverse)
library(rio)
library(table1)
library(tictoc)

# for IRT models
library(mirt)
library(psych)

# for model outputs
library(modelsummary)
library(flextable)


#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#
#### Custom Functions ####
#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#

# Custom Functions
range01 <- function(x, ...){(x - min(x, na.rm = TRUE)) / (max(x, na.rm = TRUE) - min(x, na.rm = TRUE))}  ## formula adjusted to ignore NAs


#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#
#### Set Seed ####
#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#

set.seed(515)

#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#
#### Load Data ####
#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#

### CAN ###

load(here("1 Original Data", "Original Data (Canada).RData"))

### JAP ###

load(here("1 Original Data", "Original Data (Japan).RData"))


#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#
#### Measurement Analysis ####
#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#


#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#
#### ~ Factor Analysis ~ ####
#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#


#### ~~ Japan ~~ ####

controlJPN <- df_jap %>%
  select(ResponseId, contains("q10")) %>%
  mutate(across(contains("q10"), as.numeric)) %>%
  mutate(across(q10_1, ~recode(.,
                               '4' = 1,
                               '3' = 2,
                               '2' = 3,
                               '1' = 4
  )))  %>%
  mutate(across(q10_5, ~recode(.,
                               '4' = 1,
                               '3' = 2,
                               '2' = 3,
                               '1' = 4
  )))  %>%
  mutate(across(q10_7, ~recode(.,
                               '4' = 1,
                               '3' = 2,
                               '2' = 3,
                               '1' = 4
  )))


dreadJPN <- df_jap %>%
  select(ResponseId, contains("q9")) %>%
  mutate(across(contains("q9"), as.numeric)) %>%
  mutate(across(q9_1:q9_7, ~recode(.,
                                   '4' = 1,
                                   '3' = 2,
                                   '2' = 3,
                                   '1' = 4
  )))


#efaJPN <- as.data.frame(cbind(controlJPN, dreadJPN))
efaJPN <- left_join(controlJPN, dreadJPN)

head(efaJPN)


efaJPN <- efaJPN[complete.cases(efaJPN), ] ## Remove NAs
psych::scree(efaJPN[,-1], pc = FALSE)


png(file = here("4 Figures", "SI_Scree_Japan.png"),
    #type = "cairo",
    units="in", 
    width=10, 
    height=8, 
    pointsize=12, 
    res=96)

psych::scree(efaJPN[,-1], pc = FALSE)

dev.off()


## According to the Kaiser criterion, Nfacs = 2
Nfacs <- 2
fit_efaJPN <- factanal(efaJPN[,-1], Nfacs, rotation="promax")
print(fit_efaJPN, digits=2, cutoff=0.3, sort=TRUE)
print(fit_efaJPN, digits=2, cutoff=0.3, sort=FALSE)

# remove control items with a factor loading <0.3 for both factors
#controlJPNNew <- df_jap %>%
#  select(c(q10_2,q10_3,q10_4,q10_6)) %>%
#  mutate(across(where(is.character), as.numeric))


#efaJPNnew <- as.data.frame(cbind(controlJPNNew, dreadJPN))

#head(efaJPNnew)


#efaJPNnew <- efaJPNnew[complete.cases(efaJPNnew), ] ## Remove NAs
#psych::scree(efaJPNnew, pc = FALSE)

## According to the Kaiser criterion, Nfacs = 2
#Nfacs <- 2
#fit_efaJPNnew <- factanal(efaJPNnew, Nfacs, rotation="promax")
#print(fit_efaJPNnew, digits=2, cutoff=0.3, sort=TRUE)
#print(fit_efaJPNnew, digits=2, cutoff=0.3, sort=FALSE)

# scree plot and factor analysis are largely the same after removing the control items


#### ~~ Canada ~~ ####

controlCDN <- df_can %>%
  select(ResponseId, contains("q11")) %>%
  mutate(across(contains("q11"), as.numeric)) %>%
  mutate(across(q11_1:q11_7, ~ replace(., . %in% c(-99), NA)))  %>%
  mutate(across(q11_1, ~recode(.,
                               '4' = 1,
                               '3' = 2,
                               '2' = 3,
                               '1' = 4
  )))  %>%
  mutate(across(q11_5, ~recode(.,
                               '4' = 1,
                               '3' = 2,
                               '2' = 3,
                               '1' = 4
  )))  %>%
  mutate(across(q11_7, ~recode(.,
                               '4' = 1,
                               '3' = 2,
                               '2' = 3,
                               '1' = 4
  )))

dreadCDN <- df_can %>%
  select(ResponseId, contains("q10")) %>%
  mutate(across(contains("q10"), as.numeric)) %>%
  mutate(across(q10_1:q10_7, ~ replace(., . %in% c(-99), NA)))  %>%
  mutate(across(q10_1:q10_7, ~recode(.,
                                     '5' = 1,
                                     '4' = 2,
                                     '3' = 3,
                                     '2' = 4
  )))

#efaCDN <- as.data.frame(cbind(controlCDN, dreadCDN))
efaCDN <- left_join(controlCDN, dreadCDN)

head(efaCDN)


efaCDN <- efaCDN[complete.cases(efaCDN), ] ## Remove NAs
psych::scree(efaCDN[,-1], pc = FALSE)

png(file = here("4 Figures", "SI_Scree_Canada.png"),
    #type = "cairo",
    units="in", 
    width=10, 
    height=8, 
    pointsize=12, 
    res=96)

psych::scree(efaCDN[,-1], pc = FALSE)

dev.off()


## According to the Kaiser criterion, Nfacs = 2
Nfacs <- 2
fit_efaCDN <- factanal(efaCDN[,-1], Nfacs, rotation="promax")
print(fit_efaCDN, digits=2, cutoff=0.3, sort=TRUE)
print(fit_efaCDN, digits=2, cutoff=0.3, sort=FALSE)



#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#
#### ~ IRT ~ ####
#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#

#### ~~ Control ~~ ####

can_4item_control <- controlCDN %>%
  select(
    ResponseId,
    q2=q11_2, 
    q3=q11_3, 
    q4=q11_4, 
    q6=q11_6
    ) %>%
  mutate(country = "Canada") %>%
  drop_na()

jap_4item_control <- controlJPN %>%
  select(
    ResponseId,
    q2=q10_2, 
    q3=q10_3, 
    q4=q10_4, 
    q6=q10_6
  ) %>%
  mutate(country = "Japan") %>%
  drop_na()


# graded response model (Canada - Control)
grm_can_4item_control <- mirt(can_4item_control[,2:5], 1, itemtype = "graded", SE = TRUE)
coef(grm_can_4item_control, IRTpars = TRUE)

#plot(grm_can_4item_control, type = "trace", which.items = 1:ncol(can_4item_control[,2:5]))
plot(grm_can_4item_control, type = "infoSE", which.items = 1:ncol(can_4item_control[,2:5]))
plot(grm_can_4item_control, type = 'infotrace')

summary(grm_can_4item_control)

## Calculate reliability of estimated
## theta scores. 

grm_can_4item_control_scores <- fscores(grm_can_4item_control, full.scores.SE = TRUE)

empirical_rxx(grm_can_4item_control_scores)
#0.75

# graded response model (Japan - Control)
grm_jap_4item_control <- mirt(jap_4item_control[,2:5], 1, itemtype = "graded", SE = TRUE)
coef(grm_jap_4item_control, IRTpars = TRUE)

#plot(grm_jap_4item_control, type = "trace", which.items = 1:ncol(jap_4item_control[,2:5]))
plot(grm_jap_4item_control, type = "infoSE", which.items = 1:ncol(jap_4item_control[,2:5]))
plot(grm_jap_4item_control, type = 'infotrace')

summary(grm_jap_4item_control)

## Calculate reliability of estimated
## theta scores. 

grm_jap_4item_control_scores <- fscores(grm_jap_4item_control, full.scores.SE = TRUE)

empirical_rxx(grm_jap_4item_control_scores)
#0.72


#### ~~~ Make Final Control Measures ~~~ ####

# Make data Canada
head(can_4item_control)
can_4item_control$control_srm4 <- range01(apply(can_4item_control[,2:5], 1, mean))
psych::alpha(scale(can_4item_control[,2:5]))
# alpha = 0.71
can_4item_control$control_grm4 <- range01(fscores(grm_can_4item_control, full.scores.SE = TRUE)[,1])
can_4item_control$control_bestq2 <- range01(can_4item_control$q2)


###

can_control_pairsdata <- can_4item_control %>%
  select(
    'Summated Rating Model'=control_srm4,
    'Graded Response Model'=control_grm4,
    'Single Question'=control_bestq2
  )


GGally::ggpairs(can_control_pairsdata) +
  theme_bw()

png(file = here("4 Figures", "SI_Corr_Canada_Control.png"),
    type = "cairo",
    units="in", 
    width=10, 
    height=8, 
    pointsize=12, 
    res=96)

GGally::ggpairs(can_control_pairsdata) +
  theme_bw()

dev.off()

###


can_4item_control <- can_4item_control %>%
  select(ResponseId, country, contains("control"))
head(can_4item_control)


# Make data
head(jap_4item_control)
jap_4item_control$control_srm4 <- range01(apply(jap_4item_control[,2:5], 1, mean))
psych::alpha(scale(jap_4item_control[,2:5]))
# alpha = 0.7
jap_4item_control$control_grm4 <- range01(fscores(grm_jap_4item_control, full.scores.SE = TRUE)[,1])
jap_4item_control$control_bestq2 <- range01(jap_4item_control$q2)

###

jap_control_pairsdata <- jap_4item_control %>%
  select(
    'Summated Rating Model'=control_srm4,
    'Graded Response Model'=control_grm4,
    'Single Question'=control_bestq2
  )


GGally::ggpairs(jap_control_pairsdata) +
  theme_bw()

png(file = here("4 Figures", "SI_Corr_Japan_Control.png"),
    type = "cairo",
    units="in", 
    width=10, 
    height=8, 
    pointsize=12, 
    res=96)

GGally::ggpairs(jap_control_pairsdata) +
  theme_bw()

dev.off()

###

jap_4item_control <- jap_4item_control %>%
  select(ResponseId, country, contains("control"))
head(jap_4item_control)



#### ~~ Dread ~~ ####


can_4item_dread <- dreadCDN %>%
  select(
    ResponseId,
    q7=q10_7, 
    q6=q10_6, 
    q5=q10_5, 
    q3=q10_3
  ) %>%
  mutate(country = "Canada") %>%
  drop_na()

jap_4item_dread <- dreadJPN %>%
  select(
    ResponseId,
    q7=q9_7, 
    q6=q9_6, 
    q5=q9_5, 
    q3=q9_3
  ) %>%
  mutate(country = "Japan") %>%
  drop_na()


# graded response model (Canada - Dread)
grm_can_4item_dread <- mirt(can_4item_dread[,2:5], 1, itemtype = "graded", SE = TRUE)
coef(grm_can_4item_dread, IRTpars = TRUE)

plot(grm_can_4item_dread, type = "trace", which.items = 1:ncol(can_4item_dread[,2:5]))
plot(grm_can_4item_dread, type = "infoSE", which.items = 1:ncol(can_4item_dread[,2:5]))
plot(grm_can_4item_dread, type = 'infotrace')

summary(grm_can_4item_dread)

## Calculate reliability of estimated
## theta scores. 

grm_can_4item_dread_scores <- fscores(grm_can_4item_dread, full.scores.SE = TRUE)

empirical_rxx(grm_can_4item_dread_scores)
#0.86


# graded response model (Japan - Dread)
grm_jap_4item_dread <- mirt(jap_4item_dread[,2:5], 1, itemtype = "graded", SE = TRUE)
coef(grm_jap_4item_dread, IRTpars = TRUE)

plot(grm_jap_4item_dread, type = "trace", which.items = 1:ncol(jap_4item_dread[,2:5]))
plot(grm_jap_4item_dread, type = "infoSE", which.items = 1:ncol(jap_4item_dread[,2:5]))
plot(grm_jap_4item_dread, type = 'infotrace')

summary(grm_jap_4item_dread)

## Calculate reliability of estimated
## theta scores. 

grm_jap_4item_dread_scores <- fscores(grm_jap_4item_dread, full.scores.SE = TRUE)

empirical_rxx(grm_jap_4item_dread_scores)
#0.86


#### ~~~ Make Final Dread Measures ~~~ ####


# Make data
head(can_4item_dread)
can_4item_dread$dread_srm4 <- range01(apply(can_4item_dread[,2:5], 1, mean))
psych::alpha(scale(can_4item_dread[,2:5]))
# alpha = 0.86
can_4item_dread$dread_grm4 <- range01(fscores(grm_can_4item_dread, full.scores.SE = TRUE)[,1])
can_4item_dread$dread_bestq7 <- range01(can_4item_dread$q7)

###

can_dread_pairsdata <- can_4item_dread %>%
  select(
    'Summated Rating Model'=dread_srm4,
    'Graded Response Model'=dread_grm4,
    'Single Question'=dread_bestq7
  )


GGally::ggpairs(can_dread_pairsdata) +
  theme_bw()

png(file = here("4 Figures", "SI_Corr_Canada_Dread.png"),
    type = "cairo",
    units="in", 
    width=10, 
    height=8, 
    pointsize=12, 
    res=96)

GGally::ggpairs(can_dread_pairsdata) +
  theme_bw()

dev.off()

###

can_4item_dread <- can_4item_dread %>%
  select(ResponseId, country, contains("dread"))


# Make data
head(jap_4item_dread)
jap_4item_dread$dread_srm4 <- range01(apply(jap_4item_dread[,2:5], 1, mean))
psych::alpha(scale(jap_4item_dread[,2:5]))
# alpha = 0.86
jap_4item_dread$dread_grm4 <- range01(fscores(grm_jap_4item_dread, full.scores.SE = TRUE)[,1])
jap_4item_dread$dread_bestq7 <- range01(jap_4item_dread$q7)

###

jap_dread_pairsdata <- jap_4item_dread %>%
  select(
    'Summated Rating Model'=dread_srm4,
    'Graded Response Model'=dread_grm4,
    'Single Question'=dread_bestq7
  )


GGally::ggpairs(jap_dread_pairsdata) +
  theme_bw()

png(file = here("4 Figures", "SI_Corr_Japan_Dread.png"),
    type = "cairo",
    units="in", 
    width=10, 
    height=8, 
    pointsize=12, 
    res=96)

GGally::ggpairs(jap_dread_pairsdata) +
  theme_bw()

dev.off()

###

jap_4item_dread <- jap_4item_dread %>%
  select(ResponseId, country, contains("dread"))



#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#
#### Predictive Models (CAN) ####
#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#


#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#
#### ~ Covariates ####
#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#

df_can_covars <- df_can %>%
  mutate(
    # gender
    gender = case_when(q1 == "1" ~ "Man",
                       q1 == "2" ~ "Woman",
                       q1 == "3" ~ "Woman"),
    gender = fct_relevel(gender, "Man"),
    # age
    q2 = as.numeric(q2),
    age = case_when(q2 >= 18 & q2 <= 24 ~ "18-24",
                    q2 >= 25 & q2 <= 34 ~ "25-34",
                    q2 >= 35 & q2 <= 44 ~ "35-44",
                    q2 >= 45 & q2 <= 54 ~ "45-54",
                    q2 >= 55 & q2 <= 64 ~ "55-64",
                    q2 >= 65 ~ "65+"),
    age = fct_relevel(age, "18-24", "25-34", "35-44", "45-54", "55-64", "65+"),
    # university degree
    unidegree = case_when(q6 == "9" | q6 == "10" |q6 == "11" ~ "Yes",
                          TRUE ~ "No"),
    unidegree = fct_relevel(unidegree, "No"),
    # no religion
    religion = case_when(q24 == "1" ~ "No",
                         TRUE ~ "Yes"),
    religion = fct_relevel(religion, "Yes"),
    # blue collar
    bluecollar = case_when(q27 == "2" | q27 == "3" | q27 == "4" | q27 == "5" | q27 == "9" | q27 == "10" ~ "Yes",
                           TRUE ~ "No"),
    bluecollar = fct_relevel(bluecollar, "No"),
    # trust in scientists
    trustscience = case_when(q13 == "1" ~ 4,
                             q13 == "2" ~ 3,
                             q13 == "3" ~ 2,
                             q13 == "4" ~ 1),
    trustscience = range01(trustscience),
    # tech job prospects (1 = negative, 0 = positive)
    techjob = as.numeric(q28),
    techjob = range01(techjob),
    # ideology
    ideo = as.numeric(q25),
    ideo = case_when(ideo == 12 ~ 6,
                     TRUE ~ ideo),
    ideo = range01(ideo),
    # agreeableness (higher values = less agreeable)
    agree = range01(as.numeric(q8_2)),
    # emotional stability (higher values = less stable)
    emotionalstab = range01(as.numeric(q8_4))
  )

#table(df_can_covars$q1,df_can_covars$gender)
#table(df_can_covars$q6,df_can_covars$unidegree)
#table(df_can_covars$q24,df_can_covars$religion)
#table(df_can_covars$q27,df_can_covars$bluecollar)
#table(df_can_covars$q2,df_can_covars$age)
#table(df_can_covars$q13,df_can_covars$trustscience)
#table(df_can_covars$q25, df_can_covars$ideo)
#table(df_can_covars$q28, df_can_covars$techjob)
#table(df_can$q8_2)

#table(df_can_covars$q1)
#table(df_can_covars$q13)
#table(df_can_covars$q25)
#table(df_can_covars$q21_1)
#table(df_can_covars$q21_2)
#table(df_can_covars$q21_3)
#table(df_can_covars$q21_4)


#### ~~ Scale: Conspiratorial Thinking ~~ ####

df_can_conspiracy <- df_can %>%
  select(ResponseId, contains("q21")) %>%
  mutate(dplyr::across(q21_1:q21_4, ~dplyr::recode(.,
                               '4' = "1",
                               '3' = "2",
                               '2' = "3",
                               '1' = "4"
  ))) %>%
  mutate(dplyr::across(contains("q21"), as.numeric)) %>%
  drop_na()
  


head(df_can_conspiracy)

psych::alpha(scale(df_can_conspiracy[,2:5]))
# alpha = 0.84

df_can_conspiracy$conspiracyscale <- range01(apply(df_can_conspiracy[,2:5], 1, mean))

df_can_conspiracy <- df_can_conspiracy %>%
  select(ResponseId, conspiracyscale)

# NOTE higher values = more conspiratorial


#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#
#### ~ Clean Data ~ ####
#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#

df_can_analysis <- df_can_covars %>%
  select(ResponseId, gender, age, unidegree, religion, bluecollar, trustscience, techjob, ideo, agree, emotionalstab) %>%
  left_join(df_can_conspiracy) %>%
  left_join(can_4item_control) %>%
  left_join(can_4item_dread)

head(df_can_analysis)

save(df_can_analysis, file = here("2 Clean Data", "Clean Data for Predictive Models (Canada).RData"))

#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#
#### ~ Analysis ~ ####
#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#

cortest_can <- df_can_analysis %>%
  select(dread_srm4, control_srm4) %>%
  drop_na()

cor(cortest_can$dread_srm4, cortest_can$control_srm4, method = "pearson")
#0.05189031


mod_can_control_srm4 <- lm(control_srm4 ~ gender + unidegree + religion + bluecollar + age + 
                             trustscience + techjob + ideo + conspiracyscale + agree + emotionalstab, data = df_can_analysis)
summary(mod_can_control_srm4)

mod_can_control_grm4 <- lm(control_grm4 ~ gender + unidegree + religion + bluecollar + age + 
                             trustscience + techjob + ideo + conspiracyscale + agree + emotionalstab, data = df_can_analysis)
summary(mod_can_control_grm4)

mod_can_control_bestq2 <- lm(control_bestq2 ~ gender + unidegree + religion + bluecollar + age + 
                               trustscience + techjob + ideo + conspiracyscale + agree + emotionalstab, data = df_can_analysis)
summary(mod_can_control_bestq2)



mod_can_dread_srm4 <- lm(dread_srm4 ~ gender + unidegree + religion + bluecollar + age + 
                           trustscience + techjob + ideo + conspiracyscale + agree + emotionalstab, data = df_can_analysis)
summary(mod_can_dread_srm4)

mod_can_dread_grm4 <- lm(dread_grm4 ~ gender + unidegree + religion + bluecollar + age + 
                           trustscience + techjob + ideo + conspiracyscale + agree + emotionalstab, data = df_can_analysis)
summary(mod_can_dread_grm4)

mod_can_dread_bestq7 <- lm(dread_bestq7 ~ gender + unidegree + religion + bluecollar + age + 
                             trustscience + techjob + ideo + conspiracyscale + agree + emotionalstab, data = df_can_analysis)
summary(mod_can_dread_bestq7)






#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#
#### Predictive Models (JAP) ####
#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#


#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#
#### ~ Covariates ####
#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#

df_jap_covars <- df_jap %>%
  mutate(
    # gender
    gender = case_when(q1 == "1" ~ "Man",
                       q1 == "2" ~ "Woman",
                       q1 == "3" ~ "Woman"),
    gender = fct_relevel(gender, "Man"),
    # age
    q2 = as.numeric(q2),
    q2 = 2024-q2-1,
    age = case_when(q2 >= 18 & q2 <= 24 ~ "18-24",
                    q2 >= 25 & q2 <= 34 ~ "25-34",
                    q2 >= 35 & q2 <= 44 ~ "35-44",
                    q2 >= 45 & q2 <= 54 ~ "45-54",
                    q2 >= 55 & q2 <= 64 ~ "55-64",
                    q2 >= 65 ~ "65+"),
    age = fct_relevel(age, "18-24", "25-34", "35-44", "45-54", "55-64", "65+"),
    # university degree
    unidegree = case_when(q5 == "5" ~ "Yes",
                          TRUE ~ "No"),
    unidegree = fct_relevel(unidegree, "No"),
    # no religion
    religion = case_when(q21 == "10" ~ "No",
                         TRUE ~ "Yes"),
    religion = fct_relevel(religion, "Yes"),
    # blue collar
    bluecollar = case_when(q24 == "2" | q24 == "3" | q24 == "4" | q24 == "5" | q24 == "9" | q24 == "10" ~ "Yes",
                           TRUE ~ "No"),
    bluecollar = fct_relevel(bluecollar, "No"),
    # trust in scientists
    trustscience = case_when(q12 == "1" ~ 4,
                             q12 == "2" ~ 3,
                             q12 == "3" ~ 2,
                             q12 == "4" ~ 1),
    trustscience = range01(trustscience),
    # tech job prospects (1 = negative, 0 = positive)
    techjob = as.numeric(q25),
    techjob = range01(techjob),
    # ideology
    ideo = as.numeric(q22),
    ideo = case_when(ideo == 12 ~ 6,
                     TRUE ~ ideo),
    ideo = replace_na(ideo, 6),
    ideo = range01(ideo),
    # agreeableness (higher values = less agreeable)
    agree = range01(as.numeric(q7_2)),
    # emotional stability (higher values = less stable)
    emotionalstab = range01(as.numeric(q7_4))
    
  )

#table(df_jap_covars$q1,df_jap_covars$gender)
#table(df_jap_covars$q6,df_jap_covars$unidegree)
#table(df_jap_covars$q24,df_jap_covars$religion)
#table(df_jap_covars$q27,df_jap_covars$bluecollar)
#table(df_jap_covars$q2,df_jap_covars$age)
#table(df_jap_covars$q13,df_jap_covars$trustscience)
#table(df_jap_covars$q25, df_jap_covars$ideo)
#table(df_jap_covars$q28, df_jap_covars$techjob)

#### ~~ Scale: Conspiratorial Thinking ~~ ####

df_jap_conspiracy <- df_jap %>%
  select(ResponseId, contains("q18")) %>%
  mutate(dplyr::across(q18_1:q18_4, ~dplyr::recode(.,
                                                   '4' = "1",
                                                   '3' = "2",
                                                   '2' = "3",
                                                   '1' = "4"
  ))) %>%
  mutate(dplyr::across(contains("q18"), as.numeric)) %>%
  drop_na()



head(df_jap_conspiracy)

psych::alpha(scale(df_jap_conspiracy[,2:5]))
# alpha = 0.83

df_jap_conspiracy$conspiracyscale <- range01(apply(df_jap_conspiracy[,2:5], 1, mean))

df_jap_conspiracy <- df_jap_conspiracy %>%
  select(ResponseId, conspiracyscale)

# NOTE higher values = more conspiratorial


#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#
#### ~ Clean Data ~ ####
#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#

df_jap_analysis <- df_jap_covars %>%
  select(ResponseId, gender, age, unidegree, religion, bluecollar, trustscience, techjob, ideo, agree, emotionalstab) %>%
  left_join(df_jap_conspiracy) %>%
  left_join(jap_4item_control) %>%
  left_join(jap_4item_dread)

head(df_jap_analysis)


save(df_jap_analysis, file = here("2 Clean Data", "Clean Data for Predictive Models (Japan).RData"))


#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#
#### ~ Analysis ~ ####
#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#

cortest_jap <- df_jap_analysis %>%
  select(dread_srm4, control_srm4) %>%
  drop_na()

cor(cortest_jap$dread_srm4, cortest_jap$control_srm4, method = "pearson")
#0.0099

mod_jap_control_srm4 <- lm(control_srm4 ~ gender + unidegree + religion + bluecollar + age + 
                             trustscience + techjob + ideo + conspiracyscale + agree + emotionalstab, data = df_jap_analysis)
summary(mod_jap_control_srm4)

mod_jap_control_grm4 <- lm(control_grm4 ~ gender + unidegree + religion + bluecollar + age + 
                             trustscience + techjob + ideo + conspiracyscale + agree + emotionalstab, data = df_jap_analysis)
summary(mod_jap_control_grm4)

mod_jap_control_bestq2 <- lm(control_bestq2 ~ gender + unidegree + religion + bluecollar + age + 
                                       trustscience + techjob + ideo + conspiracyscale + agree + emotionalstab, data = df_jap_analysis)
summary(mod_jap_control_bestq2)



mod_jap_dread_srm4 <- lm(dread_srm4 ~ gender + unidegree + religion + bluecollar + age + 
                           trustscience + techjob + ideo + conspiracyscale + agree + emotionalstab, data = df_jap_analysis)
summary(mod_jap_dread_srm4)

mod_jap_dread_grm4 <- lm(dread_grm4 ~ gender + unidegree + religion + bluecollar + age + 
                           trustscience + techjob + ideo + conspiracyscale + agree + emotionalstab, data = df_jap_analysis)
summary(mod_jap_dread_grm4)

mod_jap_dread_bestq7 <- lm(dread_bestq7 ~ gender + unidegree + religion + bluecollar + age + 
                             trustscience + techjob + ideo + conspiracyscale + agree + emotionalstab, data = df_jap_analysis)
summary(mod_jap_dread_bestq7)




#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#
#### Predictive Model Figures ####
#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#

detach("package:psych", unload=TRUE)
library(tidyverse)

#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#
#### ~ Controllability ~ ####
#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#

plotdat_can_control_srm4 <- broom::tidy(mod_can_control_srm4, conf.int = TRUE) %>%
  filter(term != "(Intercept)") %>%
  mutate(statsig = case_when(conf.low > 0 | conf.high < 0 ~ "Yes",
                             TRUE ~ "No")) %>%
  mutate(termpretty = case_when(term == "genderWoman" ~ "Gender (Woman)",
                                term == "unidegreeYes" ~ "Uni. Degree (Yes)",
                                term == "religionNo" ~ "Religion (No)",
                                term == "bluecollarYes" ~ "Blue Collar (Yes)",
                                term == "age25-34" ~ "Age (25-34)",
                                term == "age35-44" ~ "Age (35-44)",
                                term == "age45-54" ~ "Age (45-54)",
                                term == "age55-64" ~ "Age (55-64)",
                                term == "age65+" ~ "Age (65+)",
                                term == "trustscience" ~ "Trust in Scientists (0-1)",
                                term == "techjob" ~ "Tech Job Prospects (0-1)",
                                term == "ideo" ~ "Left-Right Ideology (0-1)",
                                term == "conspiracyscale" ~ "Conspiracy Thinking (0-1)",
                                term == "agree" ~ "Agreeableness (0-1)",
                                term == "emotionalstab" ~ "Emotional Stability (0-1)")) %>%
  mutate(Outcome = "4-Item Scale") %>%
  mutate(Country = "Canada") %>%
  mutate(Category = case_when(term == "conspiracyscale" | 
                                term == "ideo" | 
                                term == "techjob" |
                                term == "agree" |
                                term == "emotionalstab" |
                                term == "trustscience" ~ "Predispositions",
                              TRUE ~ "Sociodemographics"))

head(plotdat_can_control_srm4)

plotdat_can_control_1q <- broom::tidy(mod_can_control_bestq2, conf.int = TRUE) %>%
  filter(term != "(Intercept)") %>%
  mutate(statsig = case_when(conf.low > 0 | conf.high < 0 ~ "Yes",
                             TRUE ~ "No")) %>%
  mutate(termpretty = case_when(term == "genderWoman" ~ "Gender (Woman)",
                                term == "unidegreeYes" ~ "Uni. Degree (Yes)",
                                term == "religionNo" ~ "Religion (No)",
                                term == "bluecollarYes" ~ "Blue Collar (Yes)",
                                term == "age25-34" ~ "Age (25-34)",
                                term == "age35-44" ~ "Age (35-44)",
                                term == "age45-54" ~ "Age (45-54)",
                                term == "age55-64" ~ "Age (55-64)",
                                term == "age65+" ~ "Age (65+)",
                                term == "trustscience" ~ "Trust in Scientists (0-1)",
                                term == "techjob" ~ "Tech Job Prospects (0-1)",
                                term == "ideo" ~ "Left-Right Ideology (0-1)",
                                term == "conspiracyscale" ~ "Conspiracy Thinking (0-1)",
                                term == "agree" ~ "Agreeableness (0-1)",
                                term == "emotionalstab" ~ "Emotional Stability (0-1)")) %>%
  mutate(Outcome = "Single Question") %>%
  mutate(Country = "Canada") %>%
  mutate(Category = case_when(term == "conspiracyscale" | 
                                term == "ideo" | 
                                term == "techjob" |
                                term == "agree" |
                                term == "emotionalstab" |
                                term == "trustscience" ~ "Predispositions",
                              TRUE ~ "Sociodemographics"))

head(plotdat_can_control_1q)

plotdat_jap_control_srm4 <- broom::tidy(mod_jap_control_srm4, conf.int = TRUE) %>%
  filter(term != "(Intercept)") %>%
  mutate(statsig = case_when(conf.low > 0 | conf.high < 0 ~ "Yes",
                             TRUE ~ "No")) %>%
  mutate(termpretty = case_when(term == "genderWoman" ~ "Gender (Woman)",
                                term == "unidegreeYes" ~ "Uni. Degree (Yes)",
                                term == "religionNo" ~ "Religion (No)",
                                term == "bluecollarYes" ~ "Blue Collar (Yes)",
                                term == "age25-34" ~ "Age (25-34)",
                                term == "age35-44" ~ "Age (35-44)",
                                term == "age45-54" ~ "Age (45-54)",
                                term == "age55-64" ~ "Age (55-64)",
                                term == "age65+" ~ "Age (65+)",
                                term == "trustscience" ~ "Trust in Scientists (0-1)",
                                term == "techjob" ~ "Tech Job Prospects (0-1)",
                                term == "ideo" ~ "Left-Right Ideology (0-1)",
                                term == "conspiracyscale" ~ "Conspiracy Thinking (0-1)",
                                term == "agree" ~ "Agreeableness (0-1)",
                                term == "emotionalstab" ~ "Emotional Stability (0-1)")) %>%
  mutate(Outcome = "4-Item Scale") %>%
  mutate(Country = "Japan") %>%
  mutate(Category = case_when(term == "conspiracyscale" | 
                                term == "ideo" | 
                                term == "techjob" |
                                term == "agree" |
                                term == "emotionalstab" |
                                term == "trustscience" ~ "Predispositions",
                              TRUE ~ "Sociodemographics"))

head(plotdat_jap_control_srm4)

plotdat_jap_control_1q <- broom::tidy(mod_jap_control_bestq2, conf.int = TRUE) %>%
  filter(term != "(Intercept)") %>%
  mutate(statsig = case_when(conf.low > 0 | conf.high < 0 ~ "Yes",
                             TRUE ~ "No")) %>%
  mutate(termpretty = case_when(term == "genderWoman" ~ "Gender (Woman)",
                                term == "unidegreeYes" ~ "Uni. Degree (Yes)",
                                term == "religionNo" ~ "Religion (No)",
                                term == "bluecollarYes" ~ "Blue Collar (Yes)",
                                term == "age25-34" ~ "Age (25-34)",
                                term == "age35-44" ~ "Age (35-44)",
                                term == "age45-54" ~ "Age (45-54)",
                                term == "age55-64" ~ "Age (55-64)",
                                term == "age65+" ~ "Age (65+)",
                                term == "trustscience" ~ "Trust in Scientists (0-1)",
                                term == "techjob" ~ "Tech Job Prospects (0-1)",
                                term == "ideo" ~ "Left-Right Ideology (0-1)",
                                term == "conspiracyscale" ~ "Conspiracy Thinking (0-1)",
                                term == "agree" ~ "Agreeableness (0-1)",
                                term == "emotionalstab" ~ "Emotional Stability (0-1)")) %>%
  mutate(Outcome = "Single Question") %>%
  mutate(Country = "Japan") %>%
  mutate(Category = case_when(term == "conspiracyscale" | 
                                term == "ideo" | 
                                term == "techjob" |
                                term == "agree" |
                                term == "emotionalstab" |
                                term == "trustscience" ~ "Predispositions",
                              TRUE ~ "Sociodemographics"))

head(plotdat_jap_control_1q)


plotdat_combo_control <- rbind(plotdat_can_control_srm4,
                               plotdat_can_control_1q,
                               plotdat_jap_control_srm4,
                               plotdat_jap_control_1q) %>%
  mutate(Category = fct_relevel(Category, "Sociodemographics")) %>%
  mutate(Outcome = fct_relevel(Outcome, "Single Question"))


plot_combo_control <- ggplot(plotdat_combo_control, aes(x = estimate, y = termpretty, xmin = conf.low, xmax = conf.high)) +
  geom_vline(xintercept = 0, color = "indianred", lty = "dashed") +
  geom_pointrange(aes(alpha = statsig, shape = Outcome), linewidth = 1.25, position=position_dodge(width=0.55)) +
  scale_alpha_manual(values = c(0.2, 1), guide=F) +
  scale_shape_manual(values = c(19,0)) +
  theme_bw() +
  theme(
    #axis.title.x = element_text(hjust = 0),
    axis.title.y = element_text(hjust = 1),
    #panel.grid.major.x = element_blank(),
    #panel.grid.minor.x = element_blank(),
    strip.text = element_text(face = "bold",
                              size = rel(0.75)),
    strip.background = element_rect(fill = "grey90", color = NA),
    legend.position = "bottom"
  )+
  guides(shape = guide_legend(reverse=T)) +
  facet_grid(Category~Country, scales = "free_y") +
  labs(
    x = "Estimate",
    y = ""
  )
  

plot_combo_control


png(file = here("4 Figures", "Main_Control.png"),
    type = "cairo",
    units="in", 
    width=10, 
    height=8, 
    pointsize=12, 
    res=96)

plot_combo_control

dev.off()



#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#
#### ~ Dread ~ ####
#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#

plotdat_can_dread_srm4 <- broom::tidy(mod_can_dread_srm4, conf.int = TRUE) %>%
  filter(term != "(Intercept)") %>%
  mutate(statsig = case_when(conf.low > 0 | conf.high < 0 ~ "Yes",
                             TRUE ~ "No")) %>%
  mutate(termpretty = case_when(term == "genderWoman" ~ "Gender (Woman)",
                                term == "unidegreeYes" ~ "Uni. Degree (Yes)",
                                term == "religionNo" ~ "Religious (No)",
                                term == "bluecollarYes" ~ "Blue Collar (Yes)",
                                term == "age25-34" ~ "Age (25-34)",
                                term == "age35-44" ~ "Age (35-44)",
                                term == "age45-54" ~ "Age (45-54)",
                                term == "age55-64" ~ "Age (55-64)",
                                term == "age65+" ~ "Age (65+)",
                                term == "trustscience" ~ "Trust in Scientists (0-1)",
                                term == "techjob" ~ "Tech Job Prospects (0-1)",
                                term == "ideo" ~ "Left-Right Ideology (0-1)",
                                term == "conspiracyscale" ~ "Conspiracy Thinking (0-1)",
                                term == "agree" ~ "Agreeableness (0-1)",
                                term == "emotionalstab" ~ "Emotional Stability (0-1)")) %>%
  mutate(Outcome = "4-Item Scale") %>%
  mutate(Country = "Canada") %>%
  mutate(Category = case_when(term == "conspiracyscale" | 
                                term == "ideo" | 
                                term == "techjob" |
                                term == "agree" |
                                term == "emotionalstab" |
                                term == "trustscience" ~ "Predispositions",
                              TRUE ~ "Sociodemographics"))

head(plotdat_can_dread_srm4)

plotdat_can_dread_1q <- broom::tidy(mod_can_dread_bestq7, conf.int = TRUE) %>%
  filter(term != "(Intercept)") %>%
  mutate(statsig = case_when(conf.low > 0 | conf.high < 0 ~ "Yes",
                             TRUE ~ "No")) %>%
  mutate(termpretty = case_when(term == "genderWoman" ~ "Gender (Woman)",
                                term == "unidegreeYes" ~ "Uni. Degree (Yes)",
                                term == "religionNo" ~ "Religious (No)",
                                term == "bluecollarYes" ~ "Blue Collar (Yes)",
                                term == "age25-34" ~ "Age (25-34)",
                                term == "age35-44" ~ "Age (35-44)",
                                term == "age45-54" ~ "Age (45-54)",
                                term == "age55-64" ~ "Age (55-64)",
                                term == "age65+" ~ "Age (65+)",
                                term == "trustscience" ~ "Trust in Scientists (0-1)",
                                term == "techjob" ~ "Tech Job Prospects (0-1)",
                                term == "ideo" ~ "Left-Right Ideology (0-1)",
                                term == "conspiracyscale" ~ "Conspiracy Thinking (0-1)",
                                term == "agree" ~ "Agreeableness (0-1)",
                                term == "emotionalstab" ~ "Emotional Stability (0-1)")) %>%
  mutate(Outcome = "Single Question") %>%
  mutate(Country = "Canada") %>%
  mutate(Category = case_when(term == "conspiracyscale" | 
                                term == "ideo" | 
                                term == "techjob" |
                                term == "agree" |
                                term == "emotionalstab" |
                                term == "trustscience" ~ "Predispositions",
                              TRUE ~ "Sociodemographics"))

head(plotdat_can_dread_1q)

plotdat_jap_dread_srm4 <- broom::tidy(mod_jap_dread_srm4, conf.int = TRUE) %>%
  filter(term != "(Intercept)") %>%
  mutate(statsig = case_when(conf.low > 0 | conf.high < 0 ~ "Yes",
                             TRUE ~ "No")) %>%
  mutate(termpretty = case_when(term == "genderWoman" ~ "Gender (Woman)",
                                term == "unidegreeYes" ~ "Uni. Degree (Yes)",
                                term == "religionNo" ~ "Religious (No)",
                                term == "bluecollarYes" ~ "Blue Collar (Yes)",
                                term == "age25-34" ~ "Age (25-34)",
                                term == "age35-44" ~ "Age (35-44)",
                                term == "age45-54" ~ "Age (45-54)",
                                term == "age55-64" ~ "Age (55-64)",
                                term == "age65+" ~ "Age (65+)",
                                term == "trustscience" ~ "Trust in Scientists (0-1)",
                                term == "techjob" ~ "Tech Job Prospects (0-1)",
                                term == "ideo" ~ "Left-Right Ideology (0-1)",
                                term == "conspiracyscale" ~ "Conspiracy Thinking (0-1)",
                                term == "agree" ~ "Agreeableness (0-1)",
                                term == "emotionalstab" ~ "Emotional Stability (0-1)")) %>%
  mutate(Outcome = "4-Item Scale") %>%
  mutate(Country = "Japan") %>%
  mutate(Category = case_when(term == "conspiracyscale" | 
                                term == "ideo" | 
                                term == "techjob" |
                                term == "agree" |
                                term == "emotionalstab" |
                                term == "trustscience" ~ "Predispositions",
                              TRUE ~ "Sociodemographics"))

head(plotdat_jap_dread_srm4)

plotdat_jap_dread_1q <- broom::tidy(mod_jap_dread_bestq7, conf.int = TRUE) %>%
  filter(term != "(Intercept)") %>%
  mutate(statsig = case_when(conf.low > 0 | conf.high < 0 ~ "Yes",
                             TRUE ~ "No")) %>%
  mutate(termpretty = case_when(term == "genderWoman" ~ "Gender (Woman)",
                                term == "unidegreeYes" ~ "Uni. Degree (Yes)",
                                term == "religionNo" ~ "Religious (No)",
                                term == "bluecollarYes" ~ "Blue Collar (Yes)",
                                term == "age25-34" ~ "Age (25-34)",
                                term == "age35-44" ~ "Age (35-44)",
                                term == "age45-54" ~ "Age (45-54)",
                                term == "age55-64" ~ "Age (55-64)",
                                term == "age65+" ~ "Age (65+)",
                                term == "trustscience" ~ "Trust in Scientists (0-1)",
                                term == "techjob" ~ "Tech Job Prospects (0-1)",
                                term == "ideo" ~ "Left-Right Ideology (0-1)",
                                term == "conspiracyscale" ~ "Conspiracy Thinking (0-1)",
                                term == "agree" ~ "Agreeableness (0-1)",
                                term == "emotionalstab" ~ "Emotional Stability (0-1)")) %>%
  mutate(Outcome = "Single Question") %>%
  mutate(Country = "Japan") %>%
  mutate(Category = case_when(term == "conspiracyscale" | 
                                term == "ideo" | 
                                term == "techjob" |
                                term == "agree" |
                                term == "emotionalstab" |
                                term == "trustscience" ~ "Predispositions",
                              TRUE ~ "Sociodemographics"))

head(plotdat_jap_dread_1q)


plotdat_combo_dread <- rbind(plotdat_can_dread_srm4,
                             plotdat_can_dread_1q,
                             plotdat_jap_dread_srm4,
                             plotdat_jap_dread_1q) %>%
  mutate(Category = fct_relevel(Category, "Sociodemographics")) %>%
  mutate(Outcome = fct_relevel(Outcome, "Single Question"))


plot_combo_dread <- ggplot(plotdat_combo_dread, aes(x = estimate, y = termpretty, xmin = conf.low, xmax = conf.high)) +
  geom_vline(xintercept = 0, color = "indianred", lty = "dashed") +
  geom_pointrange(aes(alpha = statsig, shape = Outcome), linewidth = 1.25, position=position_dodge(width=0.55)) +
  scale_alpha_manual(values = c(0.2, 1), guide=F) +
  scale_shape_manual(values = c(19,0)) +
  theme_bw() +
  theme(
    #axis.title.x = element_text(hjust = 0),
    axis.title.y = element_text(hjust = 1),
    #panel.grid.major.x = element_blank(),
    #panel.grid.minor.x = element_blank(),
    strip.text = element_text(face = "bold",
                              size = rel(0.75)),
    strip.background = element_rect(fill = "grey90", color = NA),
    legend.position = "bottom"
  )+
  guides(shape = guide_legend(reverse=T)) +
  facet_grid(Category~Country, scales = "free_y") +
  labs(
    x = "Estimate",
    y = ""
  )


plot_combo_dread


png(file = here("4 Figures", "Main_Dread.png"),
    type = "cairo",
    units="in", 
    width=10, 
    height=8, 
    pointsize=12, 
    res=96)

plot_combo_dread

dev.off()





#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#
#### Predictive Model Tables ####
#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#

modelsummary(models = list(mod_can_control_srm4, mod_can_control_grm4, mod_can_control_bestq2,
                           mod_jap_control_srm4, mod_jap_control_grm4, mod_jap_control_bestq2),
             stars = TRUE)

cm <- c( 
  'unidegreeYes' = "Uni Degree (Yes)",
  'religionNo' = "Religious (No)",
  'genderWoman' = "Gender (Woman)",
  'bluecollarYes' = "Blue Collar (Yes)",
  'age65+' = "Age (65+)",
  'age55-64' = "Age (55-64)",
  'age45-54' = "Age (45-54)",
  'age35-44' = "Age (35-44)",
  'age25-34' = "Age (25-34)",
  'trustscience' = "Trust in Scientists (0-1)",
  'techjob' = "Tech Job Prospects (0-1)",
  'ideo' = "Left_Right Ideology (0-1)",  
  'emotionalstab' = "Emotional Stability (0-1)",
  'conspiracyscale' = "Conspiracy Thinking (0-1)",
  'agree' = "Agreeableness (0-1)",

  '(Intercept)' = 'Constant'
)


rows <- tribble(
  ~term, ~term, ~term, ~term, ~term, ~term, ~term,
  'Sample','Canada','Canada','Canada','Japan','Japan','Japan',
  'DV', '4-Item SRM', '4-Item GRM', 'Single Question','4-Item SRM', '4-Item GRM', 'Single Question'
)

attr(rows, 'position') <- c(34,35)


modelsummary(models = list(mod_can_control_srm4, mod_can_control_grm4, mod_can_control_bestq2,
                           mod_jap_control_srm4, mod_jap_control_grm4, mod_jap_control_bestq2),
             #output = "latex",
             coef_map = cm, 
             stars = TRUE,
             add_rows = rows,
             gof_omit = 'F|IC|Log|Adj|edf') 


SI_tab_main_control <- modelsummary(models = list(mod_can_control_srm4, mod_can_control_grm4, mod_can_control_bestq2,
                                                  mod_jap_control_srm4, mod_jap_control_grm4, mod_jap_control_bestq2),
                            output = "latex",
                            coef_map = cm, 
                            stars = TRUE,
                            add_rows = rows,
                            gof_omit = 'F|IC|Log|Adj|edf') 

kableExtra::save_kable(SI_tab_main_control, file = here("5 Tables", "SITable_MainResultsTableControl.tex"))

modelsummary(models = list(mod_can_control_srm4, mod_can_control_grm4, mod_can_control_bestq2,
                           mod_jap_control_srm4, mod_jap_control_grm4, mod_jap_control_bestq2),
            output = here("5 Tables", "SITable_MainResultsTableControl.docx"),
            coef_map = cm, 
            stars = TRUE,
            add_rows = rows,
            gof_omit = 'F|IC|Log|Adj|edf') 







modelsummary(models = list(mod_can_dread_srm4, mod_can_dread_grm4, mod_can_dread_bestq7,
                           mod_jap_dread_srm4, mod_jap_dread_grm4, mod_jap_dread_bestq7),
             stars = TRUE)

cm <- c( 
  'unidegreeYes' = "Uni Degree (Yes)",
  'religionNo' = "Religious (No)",
  'genderWoman' = "Gender (Woman)",
  'bluecollarYes' = "Blue Collar (Yes)",
  'age65+' = "Age (65+)",
  'age55-64' = "Age (55-64)",
  'age45-54' = "Age (45-54)",
  'age35-44' = "Age (35-44)",
  'age25-34' = "Age (25-34)",
  'trustscience' = "Trust in Scientists (0-1)",
  'techjob' = "Tech Job Prospects (0-1)",
  'ideo' = "Left_Right Ideology (0-1)",  
  'emotionalstab' = "Emotional Stability (0-1)",
  'conspiracyscale' = "Conspiracy Thinking (0-1)",
  'agree' = "Agreeableness (0-1)",
  
  '(Intercept)' = 'Constant'
)


rows <- tribble(
  ~term, ~term, ~term, ~term, ~term, ~term, ~term,
  'Sample','Canada','Canada','Canada','Japan','Japan','Japan',
  'DV', '4-Item SRM', '4-Item GRM', 'Single Question','4-Item SRM', '4-Item GRM', 'Single Question'
)

attr(rows, 'position') <- c(34,35)


modelsummary(models = list(mod_can_dread_srm4, mod_can_dread_grm4, mod_can_dread_bestq7,
                           mod_jap_dread_srm4, mod_jap_dread_grm4, mod_jap_dread_bestq7),
             #output = "latex",
             coef_map = cm, 
             stars = TRUE,
             add_rows = rows,
             gof_omit = 'F|IC|Log|Adj|edf') 


SI_tab_main_dread <- modelsummary(models = list(mod_can_dread_srm4, mod_can_dread_grm4, mod_can_dread_bestq7,
                                                mod_jap_dread_srm4, mod_jap_dread_grm4, mod_jap_dread_bestq7),
                                  output = "latex",
                                  coef_map = cm, 
                                  stars = TRUE,
                                  add_rows = rows,
                                  gof_omit = 'F|IC|Log|Adj|edf') 

kableExtra::save_kable(SI_tab_main_dread, file = here("5 Tables", "SITable_MainResultsTableDread.tex"))


modelsummary(models = list(mod_can_dread_srm4, mod_can_dread_grm4, mod_can_dread_bestq7,
                           mod_jap_dread_srm4, mod_jap_dread_grm4, mod_jap_dread_bestq7),
             output = here("5 Tables", "SITable_MainResultsTableDread.docx"),
             coef_map = cm, 
             stars = TRUE,
             add_rows = rows,
             gof_omit = 'F|IC|Log|Adj|edf') 



#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#
#### Descriptive Tables ####
#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#


df_can_ds <- df_can_analysis %>%
  mutate(Country = "Canada") %>%
  rename(dread_bestq=dread_bestq7,
         control_bestq=control_bestq2)
names(df_can_ds)


label(df_can_ds$gender) <- "Gender"
label(df_can_ds$age) <- "Age"
label(df_can_ds$unidegree) <- "Uni Degree"
label(df_can_ds$religion) <- "Religion"
label(df_can_ds$bluecollar) <- "Blue Collar"
label(df_can_ds$trustscience) <- "Trust in Scientists"
label(df_can_ds$techjob) <- "Tech Job Prospects"
label(df_can_ds$ideo) <- "Left-Right Ideology"
label(df_can_ds$agree) <- "Agreeableness"
label(df_can_ds$emotionalstab) <- "Emotional Stability"
label(df_can_ds$conspiracyscale) <- "Conspiracy Thinking"

table1::table1(~ unidegree + religion + gender + bluecollar + age + trustscience + techjob + ideo + emotionalstab + conspiracyscale + agree, data = df_can_ds)




df_jap_ds <- df_jap_analysis%>%
  mutate(Country = "Japan")%>%
  rename(dread_bestq=dread_bestq7,
         control_bestq=control_bestq2)
names(df_jap_ds)


label(df_jap_ds$gender) <- "Gender"
label(df_jap_ds$age) <- "Age"
label(df_jap_ds$unidegree) <- "Uni Degree"
label(df_jap_ds$religion) <- "Religion"
label(df_jap_ds$bluecollar) <- "Blue Collar"
label(df_jap_ds$trustscience) <- "Trust in Scientists"
label(df_jap_ds$techjob) <- "Tech Job Prospects"
label(df_jap_ds$ideo) <- "Left-Right Ideology"
label(df_jap_ds$agree) <- "Agreeableness"
label(df_jap_ds$emotionalstab) <- "Emotional Stability"
label(df_jap_ds$conspiracyscale) <- "Conspiracy Thinking"

table1::table1(~ unidegree + religion + gender + bluecollar + age + trustscience + techjob + ideo + emotionalstab + conspiracyscale + agree, data = df_jap_ds)



df_combo_ds <- rbind(df_can_ds, df_jap_ds)

head(df_combo_ds)

label(df_combo_ds$gender) <- "Gender"
label(df_combo_ds$age) <- "Age"
label(df_combo_ds$unidegree) <- "Uni Degree"
label(df_combo_ds$religion) <- "Religion"
label(df_combo_ds$bluecollar) <- "Blue Collar"
label(df_combo_ds$trustscience) <- "Trust in Scientists"
label(df_combo_ds$techjob) <- "Tech Job Prospects"
label(df_combo_ds$ideo) <- "Left-Right Ideology"
label(df_combo_ds$agree) <- "Agreeableness"
label(df_combo_ds$emotionalstab) <- "Emotional Stability"
label(df_combo_ds$conspiracyscale) <- "Conspiracy Thinking"

table1::table1(~ unidegree + religion + gender + bluecollar + age + trustscience + techjob + ideo + emotionalstab + conspiracyscale + agree | Country, data = df_combo_ds)


destab_combo <- table1::table1(~ unidegree + religion + gender + bluecollar + age + trustscience + techjob + ideo + emotionalstab + conspiracyscale + agree | Country, data = df_combo_ds)

destab_combo_kable <- t1kable(destab_combo, format = "latex")

kableExtra::save_kable(destab_combo_kable, file = here("5 Tables", "SITable_Descriptives.tex"))

t1flex(destab_combo) %>% 
  save_as_docx(path=here("5 Tables", "SITable_Descriptives.docx"))


label(df_combo_ds$control_srm4) <- "4-Item Summated Rating Model (Controllability)"
label(df_combo_ds$control_grm4) <- "4-Item Graded Response Model (Controllability)"
label(df_combo_ds$control_bestq) <- "Single Question (Controllability)"

label(df_combo_ds$dread_srm4) <- "4-Item Summated Rating Model (Dread)"
label(df_combo_ds$dread_grm4) <- "4-Item Graded Response Model (Dread)"
label(df_combo_ds$dread_bestq) <- "Single Question (Dread)"


table1::table1(~ control_srm4 + control_grm4 + control_bestq + dread_srm4 + dread_grm4 + dread_bestq | Country, data = df_combo_ds)

destab_combo_scale <- table1::table1(~ control_srm4 + control_grm4 + control_bestq + dread_srm4 + dread_grm4 + dread_bestq | Country, data = df_combo_ds)

destab_combo_scale_kable <- t1kable(destab_combo_scale, format = "latex")

kableExtra::save_kable(destab_combo_scale_kable, file = here("5 Tables", "SITable_Descriptives_Scales.tex"))

t1flex(destab_combo_scale) %>% 
  save_as_docx(path=here("5 Tables", "SITable_Descriptives_Scales.docx"))


#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#
#### Session Information ####
#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#

sessionInfo()

#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#
#### END ####
#--------------------------------------------------------------------------------------------------------------------------#
#--------------------------------------------------------------------------------------------------------------------------#
